#ifndef __CORRESPONDENCE_H__
#define __CORRESPONDENCE_H__

#include <vector>
#include <fstream>
#include "alltypedefs.h"


//
//  one comparison between two regions
//
class Correspondence
{
   private:

      LabelType p_labelA;
      int p_sliceA;
      int p_xA;
      int p_yA;

      LabelType p_labelB;
      int p_sliceB;
      int p_xB;
      int p_yB;

      float p_correlation;
      float p_correlationNorm;
      float p_dist;
      float p_size;


   public:

      Correspondence()
      {
      }

      Correspondence( int sliceA, LabelType labelA,
                      int sliceB, LabelType labelB,
                      int xA, int yA, int xB, int yB,
                      float correlation, float correlationNorm, 
                      float dist, float size )
      {
         p_sliceA = sliceA;
         p_labelA = labelA;
         p_sliceB = sliceB;
         p_labelB = labelB;
         p_correlation = correlation;
         p_correlationNorm = correlationNorm;
         p_dist = dist;
         p_size = size;
         p_xA = xA;
         p_yA = yA;
         p_xB = xB;
         p_yB = yB;
      }

      Correspondence( const Correspondence& C )
      {
         p_sliceA = C.p_sliceA;
         p_labelA = C.p_labelA;
         p_sliceB = C.p_sliceB;
         p_labelB = C.p_labelB;
         p_correlation = C.p_correlation;
         p_correlationNorm = C.p_correlationNorm;
         p_dist = C.p_dist;
         p_size = C.p_size;
         p_xA = C.p_xA;
         p_yA = C.p_yA;
         p_xB = C.p_xB;
         p_yB = C.p_yB;
      }

      Correspondence& operator=( const Correspondence& C )
      {
         p_sliceA = C.p_sliceA;
         p_labelA = C.p_labelA;
         p_sliceB = C.p_sliceB;
         p_labelB = C.p_labelB;
         p_correlation = C.p_correlation;
         p_correlationNorm = C.p_correlationNorm;
         p_dist = C.p_dist;
         p_size = C.p_size;
         p_xA = C.p_xA;
         p_yA = C.p_yA;
         p_xB = C.p_xB;
         p_yB = C.p_yB;

         return *this;
      }

      // Correspondence
      double compute() const
      {
         return p_correlationNorm;
      }

      // Correspondence
      double getNorm() const
      {
         return p_correlationNorm;
      }

      //
      //  -log[ C(Ri,Rj) * exp( D(Ri,Rj)/sigma^2 ) ]
      //
      // Correspondence
      double compute( float sigma, int skip ) const
      {
         // penalize skipped slices
         double alpha =  .6 ;
         double alphaval = pow( alpha, skip - 1 );
         double beta =  (double)skip ;
         //double delta = 1.0e-46; // or 3.7e-44, to be exact 
         double delta = numeric_limits<double>::min(); 
         //cout << delta << endl ;

         double cor = p_correlationNorm * alphaval;
         //cout << p_dist <<  " " << beta  << " " << sigma << endl  ;
         double dist = gaussian( p_dist, beta * sigma );
         double c = -log( (cor * dist) +delta );

         //cout << cor << " " << dist << endl ;
         //cout << "correlation: " << p_correlationNorm << endl ;
         //cout << "cor,dist: " << cor << "," << dist << endl ;

         return c;
      }

      // Correspondence
      double gaussian( double val, double sigma ) const
      {
         return exp( -(val*val) / (sigma*sigma) ); 
         //return exp( -( (val*val) / (2*sigma*sigma) ) ); 
      }

      const LabelType labelB() const { return p_labelB; };
      int sliceB() const { return p_sliceB; };
      void xyA( int& x, int &y ) const { x = p_xA, y = p_yA; };
      void xyB( int& x, int &y ) const { x = p_xB, y = p_yB; };
      float dist() const { return p_dist; };

      // Correspondence
      void writeASCII( ofstream& fileptr ) const
      {
         fileptr << p_correlation << " "
            << p_correlationNorm << " "
            << p_dist << " "
            << p_size << " " 
            << p_xB << " "
            << p_yB << " ";
      }

      // Correspondence
      void writeBinary( ofstream& fileptr ) const
      {
         fileptr.write( (char*)&p_correlation, sizeof( p_correlation ) );
         fileptr.write( (char*)&p_correlationNorm, sizeof( p_correlationNorm ) );
         fileptr.write( (char*)&p_dist, sizeof( p_dist ) );
         fileptr.write( (char*)&p_size, sizeof( p_size ) );
         fileptr.write( (char*)&p_xB, sizeof( p_xB ) );
         fileptr.write( (char*)&p_yB, sizeof( p_yB ) );
         //cout << "ka: " << p_correlation << " " << p_correlationNorm << " " 
         //   << p_dist << " " << p_size << " " 
         //   << p_xB << " " << p_yB << endl ;
      }

      // Correspondence
      void printKey()
      {
         cout << "labelA sliceA (x,y) labelB sliceB (x,y) corr corrNorm dist size: ";
         cout << endl ;
      }

      // Correspondence
      void print()
      {
         //cout << "(labelA, sliceA, labelB, sliceB, corr, corrNorm, dist, size): ";
         cout << p_labelA << " "
              << p_sliceA << " "
              << "(" << p_xA << "," << p_yA << ") " 
              << p_labelB << " "
              << p_sliceB << " "
              << "(" << p_xB << "," << p_yB << ") " 
              << p_correlation << " "
              << p_correlationNorm << " "
              << p_dist << " "
              << p_size << " = " << 
              compute( 100, p_sliceB - p_sliceA )
              << endl;
      }

      // Correspondence
      bool exits() const
      {
         if( ( p_correlation < 0 ) || ( p_correlationNorm < 0 ) )
            return false;
         return true;
      }


};

//
//  one comparison between one region and a list of regions
//
class CorrespondenceList
{
   private:

      LabelType p_label;
      int p_slice;
      int p_x, p_y;
      vector< Correspondence > p_allCorrespondences;

   public:

      CorrespondenceList()
      {
         p_label = -1;
         p_slice = -1;
         p_x = -1; p_y = -1;
      }

      CorrespondenceList( LabelType label, int slice, 
                          int x, int y ) 
      {
         p_label = label;
         p_slice = slice;
         p_x = x;
         p_y = y;
      };

      CorrespondenceList( const CorrespondenceList& cl )
      {
         p_label = cl.p_label;
         p_slice = cl.p_slice;
         p_x = cl.p_x;
         p_y = cl.p_y;
         p_allCorrespondences = cl.p_allCorrespondences;
      }

      CorrespondenceList& operator=( const CorrespondenceList& cl )
      {
         p_label = cl.p_label;
         p_slice = cl.p_slice;
         p_x = cl.p_x;
         p_y = cl.p_y;
         p_allCorrespondences = cl.p_allCorrespondences;

         return *this;
      }

      // CorrespondenceList
      void pruneCorrespondences()
      {
         vector< Correspondence >::iterator citer =
            p_allCorrespondences.begin();
         while( citer != p_allCorrespondences.end() )
         {
            //cout << citer->getNorm() <<  endl ;
            if( (*citer).getNorm() < 0 )
               citer = p_allCorrespondences.erase( citer );

            else
               citer++;
         }
      }

      // CorrespondenceList
      void xy( int& x, int& y ) const { x = p_x, y = p_y; } ;

      // CorrespondenceList
      void writeASCII( ofstream& fileptr ) const
      {
         fileptr << p_label << " " << p_x << " " << p_y << " ";
         for( unsigned int i = 0; i < p_allCorrespondences.size(); i++ )
         {
            p_allCorrespondences[i].writeASCII( fileptr );
         }
         fileptr << endl;
      }

      // CorrespondenceList
      void writeBinary( ofstream& fileptr ) const
      {
         fileptr.write( (char*)&p_label, sizeof(p_label) );
         fileptr.write( (char*)&p_x, sizeof( p_x ) );
         fileptr.write( (char*)&p_y, sizeof( p_y ) );
         //cout << "write: " << p_label << " " << p_x << " " << p_y << endl ;
         for( unsigned int i = 0; i < p_allCorrespondences.size(); i++ )
         {
            p_allCorrespondences[i].writeBinary( fileptr );
         }
      }

      // CorrespondenceList
      void readASCII( int sliceA, int sliceB,
                      vector< LabelType > allLabelsB, 
                      vector< string > parse_line, bool prune )
      {
         int j = 3;  // skip the first label and x,y coords, we've already looked at it
         for( unsigned int i = 0; i < allLabelsB.size(); i++ )
         {
            float correlation = atof(parse_line[j].c_str()); j++;
            float correlationNorm = atof(parse_line[j].c_str()); j++;
            float dist = atof(parse_line[j].c_str()); j++;
            float size = atof(parse_line[j].c_str()); j++;
            int xB = atoi( parse_line[j].c_str() ); j++;
            int yB = atoi( parse_line[j].c_str() ); j++;

            // inlcude only correspondences that exist in the data structure
            if( prune )
            {
               if( correlation > 0 )
               {
                  //cout << "pruning" << endl ;
                  Correspondence C( p_slice, p_label,
                        sliceB, allLabelsB[i],
                        p_x, p_y, xB, yB,
                        correlation, correlationNorm,
                        dist, size );

                  p_allCorrespondences.push_back( C );
               }
            }
            else
            {
               Correspondence C( p_slice, p_label,
                     sliceB, allLabelsB[i],
                     p_x, p_y, xB, yB,
                     correlation, correlationNorm,
                     dist, size );

               p_allCorrespondences.push_back( C );
            }
         }
      }

      // CorrespondenceList
      void readBinary( int sliceA, int sliceB,
                      vector< LabelType > allLabelsB, 
                      ifstream& fileptr, bool prune )
      {
         float correlation, correlationNorm, dist, size;
         int xB, yB;

         // skip the first label and x,y coords, we've already looked at it
         //LabelType ig;
         //fileptr.read( (char*)&ig, sizeof( ig ) );
         //fileptr.read( (char*)&xB, sizeof( xB ) );
         //fileptr.read( (char*)&yB, sizeof( yB ) );
         //cout << "got: " << ig << " " << xB << " " << yB << endl ;
         for( unsigned int i = 0; i < allLabelsB.size(); i++ )
         {
            fileptr.read( (char*)&correlation, sizeof( correlation ) );
            fileptr.read( (char*)&correlationNorm, sizeof( correlationNorm ) );
            fileptr.read( (char*)&dist, sizeof( dist ) );
            fileptr.read( (char*)&size, sizeof( size ) );
            fileptr.read( (char*)&xB, sizeof( xB ) );
            fileptr.read( (char*)&yB, sizeof( yB ) );

            // inlcude only correspondences that exist in the data structure
            if( prune )
            {
               if( correlation > 0 )
               {
                  //cout << "pruning" << endl ;
                  Correspondence C( p_slice, p_label,
                        sliceB, allLabelsB[i],
                        p_x, p_y, xB, yB,
                        correlation, correlationNorm,
                        dist, size );

                  p_allCorrespondences.push_back( C );
                  //C.print();
               }
            }
            // or read everything else in
            else
            {
               Correspondence C( p_slice, p_label,
                     sliceB, allLabelsB[i],
                     p_x, p_y, xB, yB,
                     correlation, correlationNorm,
                     dist, size );

               p_allCorrespondences.push_back( C );
               //C.print();
            }
         }
      }

      // CorrespondenceList
      void add( const Correspondence& C )
      {
         p_allCorrespondences.push_back( C );
      }

      // CorrespondenceList
      const Correspondence& get( int c ) const
      {
         return p_allCorrespondences[c];
      }

      // CorrespondenceList
      vector< LabelType > getAllLabels()
      {
         vector< LabelType > alllabels;
         for( unsigned int i = 0; i < p_allCorrespondences.size(); i++ )
            alllabels.push_back( p_allCorrespondences[i].labelB() );

         return alllabels;
      }

      // CorrespondenceList
      LabelType label() { return p_label; };
      int slice() { return p_slice; };

      // CorrespondenceList
      unsigned int size() const
      {
         return p_allCorrespondences.size();
      }

      // CorrespondenceList
      bool test( LabelType l, int s ) const
      {
         //cout << "l,s: " << p_label << "," << p_slice << endl ;
         //cout << "l,s: " << l << "," << s << endl ;
         if(  ( p_label == l ) && ( p_slice == s ) )
            return true;

         return false;
      }

      // CorrespondenceList
      const Correspondence* find( const LabelType search_label ) const
      {
         for( unsigned int i = 0; i < p_allCorrespondences.size(); i++ )
         {
            //cout << p_allCorrespondences[i].labelB() << " ";
            if( p_allCorrespondences[i].labelB() == search_label )
               return &(p_allCorrespondences[i]);
         }

         //cout << "Failed to find Correspondence for label " <<
         //   search_label << " on slice,label " << p_slice << ", " 
         //   << p_label << endl ;

         return NULL;
      }


};

//
//  all the comparisons between all regions between two slices 
//
class CorrespondenceMatrix
{
   private:

      int p_sliceA, p_sliceB;
      vector< LabelType > p_allSliceALabels;
      vector< LabelType > p_allSliceBLabels;
      //vector< CorrespondenceList > p_correspondenceMatrix;
      map< LabelType, CorrespondenceList > p_correspondenceMatrix;

      static void parse(const string& str, vector<string>& tokens, const string& delimiters = " ")
      {
         // Skip delimiters at beginning.
         string::size_type lastPos = str.find_first_not_of(delimiters, 0);
         // Find first "non-delimiter".
         string::size_type pos     = str.find_first_of(delimiters, lastPos);

         while (string::npos != pos || string::npos != lastPos)
         {
            // Found a token, add it to the vector.
            tokens.push_back(str.substr(lastPos, pos - lastPos));
            // Skip delimiters.  Note the "not_of"
            lastPos = str.find_first_not_of(delimiters, pos);
            // Find next "non-delimiter"
            pos = str.find_first_of(delimiters, lastPos);
         }
      }


   public:

      CorrespondenceMatrix( int sliceA, int sliceB ) 
      { 
        p_sliceA = sliceA;
        p_sliceB = sliceB;
      };

      CorrespondenceMatrix( const CorrespondenceMatrix& cm )
      {
         p_sliceA = cm.p_sliceA;
         p_sliceB = cm.p_sliceB;
         p_allSliceALabels = cm.p_allSliceALabels;
         p_allSliceBLabels = cm.p_allSliceBLabels;
         p_correspondenceMatrix = cm.p_correspondenceMatrix;
      }

      CorrespondenceMatrix& operator=( const CorrespondenceMatrix& cm )
      {
         p_sliceA = cm.p_sliceA;
         p_sliceB = cm.p_sliceB;
         p_allSliceALabels = cm.p_allSliceALabels;
         p_allSliceBLabels = cm.p_allSliceBLabels;
         p_correspondenceMatrix = cm.p_correspondenceMatrix;

         return *this;
      }

      // CorrespondenceMatrix
      void addLabels( vector< LabelType > Alabels,
                      vector< LabelType > Blabels )
      {
         p_allSliceALabels = Alabels;
         p_allSliceBLabels = Blabels;
      }

      // CorrespondenceMatrix
      void getLabels( vector< LabelType >& Alabels,
                      vector< LabelType >& Blabels )
      {
        Alabels = p_allSliceALabels;
        Blabels = p_allSliceBLabels;
      }

      // make some test correspondences
      /*
      CorrespondenceMatrix( int s1s, vector< int > s1, 
                            int s2s, vector< int > s2 ) 
      {
         for( unsigned int i = 0; i < s1.size(); i++ )
         {
            int s1_label = s1[i];
            CorrespondenceList CL( s1_label, s1s, 0, 0 );
            for( unsigned int j = 0; j < s2.size(); j++ )
            {
               int corr = rand()%10 + 1;
               Correspondence C( s1s, s1_label, 
                                 s2s, s2[j], 
                                 s1s, s1s, 
                                 s2s, s2s,
                                 0, corr, 1, 1 );
               CL.add( C );
               cout << j << ": " ; C.print(); cout << endl ;
            }
            add( CL );

            //cout << "size all: " << p_correspondenceMatrix.size() << endl ;
         }
      };
      */

      //
      //  read in the slices for the correspondence file ONLY
      //
      // CorrespondenceMatrix
      static void getSlices( const string& filename, int& A, int& B )
      {
         ifstream fileptr;
         fileptr.open( filename.c_str() );
         if( fileptr.is_open() )
         {
            //
            //  read in the first two slices
            //  TODO - must handle binary files!!!!
            // is this a binary or an ascii file?
            int test;
            fileptr.read( (char*)&test, sizeof( test ) );
            if( test == 12345 ) // read the file in as binary
            {
               //cout << "got " << test << endl ;
               //  read in the first two slices
               fileptr.read( (char*)&A, sizeof( A ) );
               fileptr.read( (char*)&B, sizeof( B ) );
            }
            else
            {
               // read the line in as ASCII
               string line;
               vector< string > parse_line;
               fileptr.seekg( 0, ios_base::beg );
               getline( fileptr, line );
               CorrespondenceMatrix::parse( line, parse_line );
               A = atoi( parse_line[0].c_str() );
               B = atoi( parse_line[1].c_str() );
            }
         }
         else
         {
            cout << "ERROR: Can't read " << filename << endl ;
            exit(1);
         }
         fileptr.close();

      }

      //
      //  read in the correspondence file
      //
      // CorrespondenceMatrix
      CorrespondenceMatrix( string filename, bool prune = true )
      {
         ifstream fileptr;
         fileptr.open( filename.c_str(), ios::in );
         if( fileptr.is_open() )
         {
            // is this a binary or an ascii file?
            int test;
            fileptr.read( (char*)&test, sizeof( test ) );
            if( test == 12345 )
            {
               //cout << "got " << test << endl ;
               // read the file in as binary
               readBinary( fileptr, prune );
            }
            else
            {
               // read the file in as ASCII
               readASCII( fileptr, prune );
            }
         }
         else
         {
            cout << "ERROR: Can't read " << filename << endl ;
         }
         fileptr.close();

      }

      // CorrespondenceMatrix
      void readBinary( ifstream& fileptr, bool prune ) 
      {
         //cout << "reading binary" << endl ;
         //
         //  read in the first two slices
         //
         fileptr.read( (char*)&p_sliceA, sizeof( p_sliceA ) );
         fileptr.read( (char*)&p_sliceB, sizeof( p_sliceB ) );
         //cout << "slices: " << p_sliceA << " " << p_sliceB << endl ;

         //
         //  read in all the number of labels for slice B
         //
         unsigned int num;
         fileptr.read( (char*)&num, sizeof( num ) );

         //cout << "blabels: " ;
         for( unsigned int i = 0; i < num; i++ )
         {
            LabelType L;
            fileptr.read( (char*)&L, sizeof( L ) );
            p_allSliceBLabels.push_back( L );
            //cout << L << " ";
         }
         //cout << endl;

         //
         // read in all the correspondence lists
         //
         while( fileptr )
         {
            LabelType labelA;
            fileptr.read( (char*)&labelA, sizeof( labelA ) );

            if( fileptr.eof() )
               break;

            int x, y;
            fileptr.read( (char*)&x, sizeof( x ) );
            fileptr.read( (char*)&y, sizeof( y ) );
            //cout << "blah " << labelA << " " << x << " " << y << endl ;

            CorrespondenceList corrList( labelA, p_sliceA, x, y );
            corrList.readBinary( p_sliceA, p_sliceB, 
                                 p_allSliceBLabels, fileptr, prune );
            //cout << "corrList size: " << corrList.size() << endl ;
            //p_correspondenceMatrix.push_back( corrList );
            p_correspondenceMatrix[ labelA ]= corrList;
            p_allSliceALabels.push_back( labelA );

         }
         //cout << "done" << endl ;
      }

      // CorrespondenceMatrix
      void readASCII( ifstream& fileptr, bool prune )
      {
         // make sure we're at the beginning of the file
         fileptr.seekg( 0, ios_base::beg );

         //
         //  read in the first two slices
         //
         string line;
         vector< string > parse_line;
         getline( fileptr, line );
         parse( line, parse_line );
         //for( unsigned int i = 0; i < parse_line.size(); i++ )
         //   cout << parse_line[i] << " ";
         //cout << endl ;
         p_sliceA = atoi( parse_line[0].c_str() );
         p_sliceB = atoi( parse_line[1].c_str() );

         //
         //  read in all the labels for slice B
         //
         getline( fileptr, line );
         //cout << line << endl;
         parse_line.clear();
         parse( line, parse_line );
         for( unsigned int i = 0; i < parse_line.size(); i++ )
         {
            p_allSliceBLabels.push_back( atol( parse_line[i].c_str() ) );
         }

         //
         // read in all the correspondences
         //
         while( getline( fileptr, line ) )
         {
            //cout << line << endl ;
            parse_line.clear();
            parse( line, parse_line );
            LabelType labelA = atol( parse_line[0].c_str() );
            CorrespondenceList corrList( 
                  /*slice label */      labelA,
                  p_sliceA,
                  /* x */               atoi( parse_line[1].c_str()),
                  /* y */               atoi( parse_line[2].c_str()) );
            corrList.readASCII( p_sliceA, p_sliceB, p_allSliceBLabels,
                                parse_line, prune );
            //p_correspondenceMatrix.push_back( corrList );
            p_correspondenceMatrix[ labelA ]= corrList;
            p_allSliceALabels.push_back( labelA );

         }
         //cout << "done" << endl ;
      }

      // CorrespondenceMatrix
      void writeBinary( const string& prefix ) const
      {
         char filename[100];
         sprintf( filename, "%s%03i-%03i.csv", prefix.c_str(), p_sliceA, p_sliceB );
         cout << "writing to " << filename << endl ;

         ofstream fileptr;
         fileptr.open( filename, ios::out | ios::binary );
         if( fileptr.is_open() )
         {
            // write out a token for the top of the file so we know it's
            // binary...
            int binary = 12345;
            fileptr.write( (char*)&binary, sizeof( binary ) );
            //cout << "binary" << binary << endl;

            //
            //  write out the first two slice indices
            //
            fileptr.write( (char*)&p_sliceA, sizeof( p_sliceA ));
            fileptr.write( (char*)&p_sliceB, sizeof( p_sliceB ));
            //cout << "blah: " << p_sliceA << " " << p_sliceB << endl ;

            //
            //  write out all the number of labels in slice B
            //
            unsigned int num = p_allSliceBLabels.size();
            fileptr.write( (char*)&num, sizeof( num ) );
            //cout << "num: " << num << endl ;

            //
            //  write out all the labels in slice B
            //
            //cout << "labels " ;
            for( unsigned int i = 0; i < p_allSliceBLabels.size(); i++ )
            {
               fileptr.write( (char*)&p_allSliceBLabels[i], 
                                 sizeof( p_allSliceBLabels[i] ) );
               //cout << p_allSliceBLabels[i] << " ";
            }
            //cout << endl ;

            //
            //  write out all the correspondences
            //
            //p_correspondenceMatrix.begin()->second.writeBinary( fileptr );
            map< LabelType, CorrespondenceList >::const_iterator cmiter;
            for( cmiter = p_correspondenceMatrix.begin(); 
                  cmiter != p_correspondenceMatrix.end(); 
                  cmiter++ )
               cmiter->second.writeBinary( fileptr );

            fileptr.close();
         }
         else
         {
            cout << "ERROR: Unable to open " << filename << endl ;
         }

      }

      // CorrespondenceMatrix
      void writeASCII( const string& prefix ) const
      {
         char filename[100];
         sprintf( filename, "%s%03i-%03i.csv", prefix.c_str(), p_sliceA, p_sliceB );
         cout << "writing to " << filename << endl ;

         ofstream fileptr;
         fileptr.open( filename, ios::out );
         if( fileptr.is_open() )
         {
            //
            //  write out the first two slice indices
            //
            fileptr << p_sliceA << " " << p_sliceB << endl ;

            //
            //  write out all the labels in slice B
            //
            for( unsigned int i = 0; i < p_allSliceBLabels.size(); i++ )
            {
               fileptr << p_allSliceBLabels[i] << " " ;
            }
            fileptr << endl ;

            //
            //  write out all the correspondences
            //
            map< LabelType, CorrespondenceList >::const_iterator cmiter;
            for( cmiter = p_correspondenceMatrix.begin(); 
                  cmiter != p_correspondenceMatrix.end(); 
                  cmiter++ )
            {
               cmiter->second.writeASCII( fileptr );
            }
            fileptr.close();
         }
         else
         {
            cout << "ERROR: Unable to open " << filename << endl ;
         }

      }

      // CorrespondenceMatrix
      void pruneCorrespondences()
      {
         //loop over each list of correspondences and remove -1
         //correspondences
         map< LabelType, CorrespondenceList >::iterator cmiter;
         for( cmiter = p_correspondenceMatrix.begin(); 
               cmiter != p_correspondenceMatrix.end(); 
               cmiter++ )
         {
            cmiter->second.pruneCorrespondences();
         }

      }

      // CorrespondenceMatrix
      void add( const LabelType l, const CorrespondenceList& CL )
      {
         p_correspondenceMatrix[ l ] = CL;
      };

      /*
      const CorrespondenceList& get( int c ) const
      {
         return p_correspondenceMatrix[c];
      }

      LabelType getLabel( int c )
      {
         return p_correspondenceMatrix[c].label();
      }
      
      int getSlice( int c )
      {
         return p_correspondenceMatrix[c].slice();
      }
      */

      // CorrespondenceMatrix
      int getSliceA() { return p_sliceA; };
      int getSliceB() { return p_sliceB; };

      // CorrespondenceMatrix
      unsigned int size() { return p_correspondenceMatrix.size(); };

      /* needs to be updated to handle maps not vectors
      CorrespondenceList get( LabelType label, int slice )
      {
         for( unsigned int i = 0; i < p_correspondenceMatrix.size(); i++ )
            if( p_correspondenceMatrix[i].test( label, slice ) )
               return p_correspondenceMatrix[i];

      }
      */

      // CorrespondenceMatrix
      const Correspondence* getCorrespondence( const LabelType& label1, int slice1,
                                               const LabelType& label2, int slice2 )
      {
         const Correspondence* C = p_correspondenceMatrix[ label1 ].find( label2 );
         if ( C != NULL )
            return C;
         /*
         for( unsigned int i = 0; i < p_correspondenceMatrix.size(); i++ )
         {
            if( p_correspondenceMatrix[i].test( label1, slice1 ) )
            {
               const Correspondence* C = p_correspondenceMatrix[i].find( label2 );
               if ( C != NULL )
                  return C;
            }
         }*/

         //cout << "Failed to find correspondence for label,slice "
         //   << label1 << ", " << slice1 << " -> " << label2 << ", " << slice2 << endl ;
              
         return NULL;
         
      }
      
};

#endif
